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ABSTRACT 

We show that the diffusion approximation breaks down for particle acceleration at 
oblique shocks with velocities typical of young supernova remnants. Higher order 
anisotropics flatten the spectral index at quasi-parallel shocks and steepen the spectral 
index at quasi-perpendicular shocks. We compare the theory with observed spectral 
indices. 
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1 INTRODUCTION 

The well-established theory of cosmic ray (CR) acceleration 
by strong non-relativistic shocks predicts a differential en- 
ergy power law spectrum n{E)dE oc E'^dE with a spectral 
index s = 2 (Axford, Leer & Skadron 1977, Krimskii 1977, 
Bell 1978, Blandford & Ostriker 1978). However, a variety 
of spectral indices is observed. Galactic CR have a steeper 
spectrum even when allowance is made for energy-dependent 
propagation losses (Hillas 2005). Steeper spectra are also 
observed in young supernova remnants (SNR) as discussed 
below. Deviations from a s = 2 spectrum are well known 
for relativistic shocks (for example: Kirk & Heavens 1989, 
Ostrowski 1991, Baring Ellison & Jones 1993, Ruffalo 1999, 
Gieseler et al 1999, Kobayakawa, Honda & Samura 2002, 
Achterberg et al 2001). Here we show that corrections to 
theory produce a spectral index that deviates significantly 
from s = 2 even at shock velocities as low as 10,000 km 
s~^ typical of young supernova remnants. The magnitude 
of the effect depends on the ratio of the velocity of the ac- 
celerating particle to the shock velocity. Hence, it may also 
apply to the acceleration of non-relativistic particles at he- 
liospheric shocks. The spectrum can either steepen or flat- 
ten depending on the angle between the shock normal and 
the large scale upstream magnetic field. Quasi-perpendicular 
shocks (Jokipii 1982, 1987) show particularly strong steep- 
ening. The effect is small for parallel shocks where changes 
in the spectral index are second order in the ratio of the 
particle velocity to the shock velocity. Non-linear hydrody- 
namic feedback of the CR pressure onto the shock structure 
can also change the spectral index (Drury 1983, Dury & 
Volk 1981, Ellison, Baring & Jones 1996) but we do not 
take account of this in the present paper. We assume that 
the fluid velocity is uniform both upstream and downstream 
of the shock. A further assumption lying at the heart of our 
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method is that the magnetic field can be divided into (i) 
a large scale field that is uniform apart from a change of 
direction and magnitude at the shock and (ii) a small scale 
field that fluctuates on scales smaller than the CR Larmor 
radius and deflects the CR trajectories through small angles 
to cumulatively scatter the CR. This ignores magnetic field 
line wandering which can also change the CR spectral index 
(Kirk, Duffy & Gallant 1996). In our calculations we further 
assume that the fiuid is always compressed by a factor of 
four at the shock, which neglects the increase in compres- 
sion at the shock due to the different equation of state of the 
relativistic component (ratio of specific heats = 4/3) and the 
effects of finite Mach number. With these assumptions we 
isolate and concentrate on the breakdown of the diffusive 
approximation and the effect of high order anisotropies in 
the CR distribution near the shock. 



2 EQUATIONS 

Cosmic rays are thought to be accelerated to high energy by 
crossing back and forth across a shock in a first order Fermi 
process, with a fractional energy gain of order Us/c at each 
crossing where Ua is the shock velocity. A small proportion 
of the CR injected at low energy proceed to cross the shock 
many times before escaping the shock environment. Statisti- 
cally this results in a power law stretching up to a few PeV in 
the galaxy and probably approaching ZeV extragalactically. 
Here we consider CR protons of mass rUp, but the analy- 
sis applies equally well to heavier ions or to electrons. The 
shape and extent of the CR spectrum is described by the 
CR distribution function f{x,p,t) in momentum p. We lo- 
cate a planar shock at a; = 0. Upstream of the shock (a; < 0) 
the ambient plasma fiows towards the shock at the shock 
velocity Us- For a strong non-relativistic shock the plasma 
velocity is reduced to Us/A in the downstream plasma. The 
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equation for the CR distribution function is 
dt ^ dx dx^^ dpx c dx^ dpx dp 

(1) 

where / is defined in the local fluid rest frame moving at 
velocity u{x) in the x direction and relativistic effects due 

to 7u = (1 — u^/c^)^^^^ / 1 have been neglected. Vx is the 
a::-componcnt of the CR velocity v, B is the local large scale 
magnetic field, and C(/) represents CR scattering by small 
scale fiuctuations in the magnetic field. In the local fiuid rest 
frame, C{f) does not change the CR energy, but merely de- 
flects its trajectory by diffusion in angle, du/dx is zero ev- 
erywhere except at the shock where it produces CR acceler- 
ation. Building on our experience of modelling laser-plasmas 
we solve the Vlasov-Fokker-Planck (VFP) equation by ex- 
pressing the CR distribution function as a sum of spherical 
harmonics (Bell et al 2006). 
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where Px = pcosO, py — psm6cos4> and pz = psinO sm(j>. 
The VFP equation can be separated into individual moment 
equations for the evolution of fi^{x,p, t): 
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where /JJ" = {l-m){l + m+ 1)//"+' - for m > 

and pf = 2l{l + The collision term C(/) on the 

right hand side reduces to an algebraic form because spher- 
ical harmonics are eigenfunctions of the Laplacian operator 
representing diffusion in 9 and 0. The co-ordinate system 
has been chosen such that By = Q and i'(x) is the collision 
frequency. The flrst three moment equations (/ = 0, 1) are 
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where ujz = ceBz/p and ojx = ceBx/p- Supernova ejecta 
initially expand at velocities up to ~ c/5, and the SNR 
shock velocity decreases to ~ c/100 during the first 1000 
years. In the limit of small u/c the / = and / = 1 equa- 
tions dominate. In the Chapman-Enskog expansion, valid 
for scalelengths much larger than the CR scattering length, 
the magnitude of //" decreases with increasing /, such that 
/™ ~ (u/c)' fa. With this expansion, /;'" can be neglected 
for /> 1, in which case the above three equations (/ = and 
/ = 1) reduce to 
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where /o and /i are real. The term udf^/dx in the first 
equation cannot be neglected as small to order 0(tt/c) since 
fo ~ (c/m)/i , but terms involving f2^ are neglected since 
they are of order {u/c)f^. /{" can be eliminated between 
the equations resulting in a diffusion equation for /q. du/dx 
is non-zero only at the shock, so the upstream diffusive so- 
lution of these equations for the CR density fo for a general 
oblique shock is 



/o (a:) = /o (0) exp(K;x) where n - 
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The corresponding CR fluxes are 

,0/ X 3Ws ,0/ N J tif \ 3Us + i0Jx)ulz fO, X 

/i (a;) = — fo (x) and /i (a;) = - — , „ fo {x) 

(6) 

An alternative tensor notation (Johnston 1960) expresses 
the CR distribution function to flrst order in u/c as / = 
./() +fi.p/p. The tensor expansion is related to the spherical 
harmonic expansion by fx = fi,fy~ — 5R{2/i}, and fz = 
Q{2fl}. We choose the y and z directions such that By — 
0. In this geometry, fx is parallel to the shock normal, fz 
is perpendicular to the shock normal and lies in the plane 
containing the magnetic field and the shock normal, and 
is perpendicular to both the shock normal and the magnetic 
field, fi, or fx, represent CR diffusion parallel to the shock 
normal. CR drift along the magnetic field is represented by 
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a combination of /i and Both real and imaginary 

parts of /i are zero everywhere for a parallel shock. 

Downstream of the shock {x > 0), the correspond- 
ing solution to the same equations is that for a uniform 
isotropic CR distribution at rest with respect to the back- 
ground plasma: fS{x) = /o"(0), /f = = where /o°(0) 
is the CR density at the shock. Matching the upstream 
and downstream CR distribution functions at the shock re- 
quires /o to be continuous at the shock. Integrating the 
(/ = l,m = 0) moment equation across the shock re- 
quires — {u/c)pdfo /dp to be continuous at the shock. 
The term -{u/c)pdfS /dp arises from the frame transfor- 
mation across the shock. Assuming a strong shock in which 
M = Us/4 downstream, the upstream flux f1 — {us/c)pdff, /dp 
must equal the downstream flux —{us/4:c)pdfo /dp, giving 
pdfo/dp = — 4/0 from equation (6). This leads to the well- 
known power law fS oc p~*, equivalent to n{E) oc E~'^ for 
CR acceleration at a strong shock. This solution is correct in 
the limit of low shock velocity for a parallel shock in which 

By = Bz = 0. 

The /o and components of the distribution function 
can be matched across a parallel shock. However, it is not 
possible to match the solution for fl across an oblique shock 
when the higher order terms (/™ for I > 1) are neglected. 
The Chapman-Enskog expansion breaks down because of 
the abrupt gradient at the shock. In the limit of small shock 
velocity, the mismatch at an oblique shock is unimportant 
and introduces only a correction of order u/c to the solution. 
Consequently, the spectrum at an oblique shock is still /o oc 
p~^ in the low velocity limit, but the corrections of order 
u/c become important at the expansion velocities of young 
SNR. The corrections are more important for oblique shocks 
than for parallel shocks for which the corrections are of order 
(u/cf . As we show below, the CR spectrum is expected to 
deviate from p~* more strongly for oblique shocks than for 
parallel shocks. 

One reason why oblique shocks deviate more strongly 
from a spectrum is that the expansion parameter is 
more correctly u/{ccos9) rather than u/c, where d is the 
angle between the magnetic field and the shock normal. CR 
diffuse away from the shock along magnetic field lines at the 
angle 9 to the shock normal. If u > ccosO, it is impossible 
for CR to escape upstream of the shock by streaming along 
field lines. Shocks with a velocity u exceeding ccos9 are 
accounted 'superluminal', even if the shock velocity is in 
fact considerably loss than the speed of light. Corrections 
to the spectral index 7, where /o (0) oc p^'' , are especially 
strong when u/ cos approaches c. 

The aim of this paper is to calculate the CR spectral 
index at oblique shocks at shock velocities relevant to young 
SNR. Integration of the first of equations (3) from far up- 
stream to far downstream gives 



7 = 



IZo (3/0 + 6/2V5 + 2k/17c) {du/dx)dx - 3w(oo)/o°(oo) 



/_°°^ (/o° + 2/275 + nf°/c) {du/dx)dx 



(7) 

where m(oo) and fo{oo) are the values of u and /q far down- 
stream. At the shock /o (0) oc p"^. To leading order in u/c, 
and for a discontinuous shock transition 



where /g (0) is the value of f" at the shock. For a shock com- 
pression of 4 (it(cx)) = Us/ 4), -y — 4 since the downstream 
CR density is uniform (/o(oo) = /o(0)) to order u/c. We 
shall show that for an oblique shock, /o (00) / /o (0), and 7 
differs from 4 to order u/c. 

The ratio /o (00)/ /o (0) is important for determining the 
spectral index 7 because /o (0) determines the rate at which 
CR are accelerated by crossing and recrossing the shock, 
while /f?(oo) determines the rate at which CR are advected 
away from the shock downstream. The spectral index is de- 
termined by the number of times CR cross the shock before 
escaping downstream. Hence 7 depends on /o (oo)//o (0), al- 
though the values of fj and /2 at the shock can also change 
the spectral index to higher order as evident in equation (7). 



3 NUMERICAL MODEL 

For all except the highest energy CR, the acceleration time 
is shorter than the time for which the shock exists, so we 
seek a time-independent solution with dfj/^/dt = for all 
I and m. We assume that the collision frequency v due to 
CR scattering by small scale magnetic field is proportional 
to the local CR Larmor frequency: 



i'{x) oc |B(a::) 
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where B is the local magnetic field, thus making the ratio of 
the CR mean free path to the Larmor radius the same both 
upstream and downstream of the shock. By definition of v 
as a frequency for small angle scattering, the coefficients //" 
decay exponcntionally at a rate [1(1 + 1) /2]i' in the absence of 
all other effects. The large-scale magnetic field is assumed to 
be uniform in the upstream plasma, Bx ~ Bq cos6, By = 0, 
Bz = Bo sin 6 and compressed by the shock to give a uniform 
downstream field, Bx = Bo cos 9, By = 0, Sz = 4i?osin^. 
As already stated, we assume throughout that the plasma 
is compressed by a factor 4 at the shock, even at high shock 
velocity where the compression is in reality increased by the 
reduced relativistic ratio of specific heats. 

Ftom the self-similar structure of equation (2) it can 
be shown that the steady state CR distribution function is 
a power law in momentum at the shock. The shock is the 
only point in space at which du/dx is non-zero, and hence 
the only point at which f(x,p) is connected to the values of 
f{x,p + Ap) and f{x,p— Ap). If the spectral index at the 
shock is 7, pdfl^ /dp can be replaced by —^fr in equation 
(2) so that it becomes an equation with differentials in x 
but not p. p only appears in combination with the magnetic 
field as cBx/p and cBz/p as the inverse of the CR Larmor 
radius by which all distances are scaled. The objective of the 
calculation is then to find the value of 7 for which the bound- 
ary conditions are //" = for all l,m far upstream where 
x —00, and for which only fo is non-zero far downstream. 
Downstream, //" — > as a; — > 00 for I > 0. 

The steady state VFP equation is non-linear because 
of the multiplication of the two unknowns / and 7 in terms 
proportional to pdfl^ /dp = —^fi". Making the replacement 
pdfr/9p = —'yfr', we iterate to a steady state solution 
by restoring the time-derivative in equation (2), fixing the 
value of fo at the shock, and improving the value of 7 as 
equation (2) is integrated forward in time until a steady 
state is reached in which CR advect at constant density 
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through a computational boundary placed far downstream 
of the shock. During the time integration, 7 is calculated 
at each timestep from equation (7). Equation (2) is solved 
in finite difference form with //" defined at the cell centres 
for even I, and defined on the cell boundaries for odd /. All 
spatial derivatives are then naturally centre-differenced. A 
stretched spatial grid is used to give fine resolution near 
the shock, where /[" changes rapidly on scales less than 
a CR Larmor radius. A coarser resolution is used at large 
distances from the shock where //" changes on the scale- 
length of the shock precursor. The background fluid veloc- 
ity u is pre-determined and changes smoothly by a factor 
of four across a thin shock transition. B is determined by 
the flux freezing condition. At the shock, wo impose the 
form u = 0.625ms — 0.375ms tanh(x/a;s). For high shock ve- 
locity Us = c/5, the lengthscale Xs for the shock transition 
is set to 0.01 times the upstream CR Larmor radius, and 
the spacing of the computational grid at the shock is set to 
Aa: = Xa/S. As shown in figure 1(a), halving Aa:: makes very 
little difference in the case requiring the highest resolution 
(a perpendicular shock with Us = c/5). Xs and Ax are set a 
factor of two larger for low shock velocities where the natural 
scalelengths are larger. For quasi-perpendicular high veloc- 
ity shocks {Us = c/5), the spherical harmonic expansion was 
extended to Z = 11. Test calculations with a larger number, 
I = 15, made negligible difference to the results. 

With our neglect of non-linear hydrodynamic effects due 
to the CR pressure in the precursor, the shock transition 
takes place over a distance very much less than the Larmor 
radius of high energy CR. For this reason our chosen scale- 
length Xs is much smaller than the Larmor radius. However 
we cannot treat the shock as a discontinuity since du/dx 
would then be infinite which is not possible in our numer- 
ical treatment. An alternative numerical method would be 
to treat the shock as a discontinuity with a matrix relat- 
ing the values of each side of the shock, = ^"i^/;"- 
A^li^ would be a full matrix in {li,l2) with the spherical 
harmonic expansion extended to very high order. A discon- 
tinuous shock introduces high order harmonics, making this 
approach unrealistic. Smoothing the shock over a small dis- 
tance Xs terminates the harmonic expansion without signif- 
icantly changing the spectral index. However, as an addi- 
tional check on the overall accuracy of the code (Code A), 
we wrote a completely separate code (Code B) which treats 
the shock jump differently, allowing a sharper transition at 
the cost of dispensing with second order terms in u/c at the 
shock. Code B solves the reduced equations 
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and all other S/" = 0. ^2 must be included because it bal- 
ances the jump in /}' across the shock. All other terms in- 
volving du/dx are smaller to order u/c but not negligible. 
Their omission substantially simplifies the numerical inte- 
gration across the shock. Some minimal smoothing of the 
shock transition on the scale of the computational grid is 



still included to allow the spherical harmonic expansion to 
be terminated. The output from the code B is compared with 
the results from the main code (code A) in figure 2(b) for 
v/iOg = 0.1 and u/c = 0.1, which we consider to be a rep- 
resentative 'standard' set of parameters. As expected, the 
comparison shows a difference in the spectral index due to 
the neglect of second order terms in u/c, but the two curves 
are sufficiently close to give confidence that the results from 
Code A can be relied upon. 



4 RESULTS 

Equation (2) was solved for a variety of shock velocities. 
Us = c/5, c/10 & c/30, and a variety of coUisionalities, 
i^/cjj = 0.03, 0.1, 0.3 & 1.0, where ujg = ceBo/p is the 
CR Larmor frequency in the upstream plasma. The spec- 
tral indices 7 are plotted in figure 1 as a function of cosO 
where 9 is the angle between the shock normal and the mag- 
netic field. At a perpendicular shock cos 6* = 0. At a parallel 
shock cos 6' = 1. In the case of a shock advancing into a 
randomly disordered magnetic field cos^ is uniformly dis- 
tributed between and 1. Also plotted on the vertical axis 
is the synchrotron spectral index a = (7 — 3)/2. 
The results can be summarised as follows: 

1) In the case of parallel shocks, the CR spectral index de- 
viates little from 7 = 4 (a = 0.5). 

2) The quasi-parallel spectrum is flattened (7 < 4, a < 0.5). 

3) The quasi-perpendicular spectrum is steepened (7 > 4, 
a > 0.5) and CR are accelerated less efficiently. 

4) Deviation from 7 = 4 is strongest for high shock velocity, 
but significant deviations are still present at velocities as low 
as 10,000 km sec"^ 

5) Deviations from 7 = 4 are greater if the angular scatter- 
ing time is much longer than the CR Larmor gyration time 

6) In the case of a shock advancing into a magnetic field 
which is disordered on scales both larger and smaller than 
a CR Larmor radius, the spectral index averages to 7 = 4 
since the mean value of 7 averaged over cosO is close to 
7 = 4 (see figure 1). 

7) The overall spectral index can only be strongly fiattened 
if the shock is predominantly quasi-parallel. Conversely, it 
can only be strongly steepened if the shock is predominantly 
quasi-perpendicular. 

The magnitude of spectral index deviation from 7 = 4 
depends strongly on the ratio of the CR scattering frequency 
V to the CR Larmor froqency u)g. The parameter v/u)g is 
poorly known both observationally and theoretically, but 
there are reasons for assuming that for SNR shocks v/ujg 
lies in the band of values we choose for figure 1. Very low 
collision frequencies {v/ujg -C 1) can be ruled out for SNR 
shocks because CR acceleration would then be too slow for 
CR acceleration to PcV energies as required by the galactic 
CR spectrum. At the opposite limit, large collision frequen- 
cies (y/ujg > 1) can probably be ruled out since v/ojg > 1 
would correspond to fields with no large scale structure. This 
is unlikely because some form of power-law magnetic field 
turbulence spectrum can be expected which includes fluc- 
tations at wavenumbers less than, as well as greater than, 
the inverse of any particular CR Larmor radius. As a guess, 
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Figure 1. Spectral indices at three sliock velocities (c/5, c/10, c/30) and four CR scattering rates {I'/uig = 0.03, 0.1, 0.3 and 1.0). 7 is 
the CR energy spectral index, a is the corresponding electron synchrotron spectral index, o = (7 — 3)/2. 9 is the shock obliquity. At a 
perpendicular shock cos6 = 0. At a parallel shock cosfl = 1. The dotted line in 1(a) next to the curve for u/ujg = 0.03 is a calculation 
with finer spatial resolution with Aa; halved as an accuracy check on the most demanding calculation. The dotted line in 1(b) next to 
the curve for u/uig = 0.1 is a comparison calculation with code B. All other curves are calculated with Code A. 



and it is little more than that at present, it seems reason- 
able to suppose that anisotropic CR drifts {fx, fy & /z) 
might decay by a factor of two in the time it takes for a 
CR to make one gyration through an angle of 27r. This cor- 
responds to f/iUg = 0.11. It is often presumed that 'Bohm 
diffusion' applies to CR scattering, but this can variously 
and ambiguously be taken to mean f/ojg = 1, I'/ujg = 1/3 
or v/ijjg = 0.125 as in the NRL Plasma Formulary (Huba 
1994). i^/ujg = 1 would seem to be too large because drift 
anisotropies then decay by a factor of ~ exp(— 27r) — 10""^ 
during one CR gyration time, which would only be possible 
if the magnetic field responsible for small scale scattering 
is much greater than the field structured on the scale of a 
CR Larmor radius. The value of the ratio i^/ujg is a major 
uncertainty in our discussion, and we take ly/ujg = 0.1 as a 
'standard' value. 

Figure 2 displays profiles of /q, /? and fl for a par- 
allel shock {9 = 0°), an oblique shock {6 = 60°), and 
a perpendicular shock {8 = 90°). Figure 3 shows the de- 
tailed profiles in the immediate environment of the shock 



for e = 60°, 72°, & 90°. In figures 2 and 3, u/ujg = 0.1 and 
Us = c/10. The oblique shock {6 — 60°) exhibits a peak in 
the CR density at the shock as previously identified in mildly 
relativistic shocks. Geissler et al (1999) convincingly argued 
for a build-up of CR immediately upstream of the shock 
because CR approaching the shock from upstream with cer- 
tain pitch angles are unable to cross the shock into the down- 
stream region. Similar structures are evident in Monte Carlo 
calculations by Ellison et al (1996). Figure 3 is a plot of CR 
density profiles at higher resolution. At ^ = 60° an excess of 
CR is seen within one CR Larmor radius both upstream and 
downstream of the shock. According to equations (7) and 
(8), as discussed above, the excess CR density at the shock 
results in enhanced CR acceleration. The excess increases 
the rate at which CR cross the shock, and hence increases 
the probability of acceleration to high energy before being 
advected away downstream of the shock. This explains the 
flattening of the CR spectrum at quasi-parallel shocks as 
seen in figure 1. 

The CR density profile at a perpendicular shock is dif- 
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Figure 2. Spatial profiles of f^, f^, 5R{/j^}, and } for a parallel shock, an oblique shock at 8 = 60° and a peperpendicular shock. 

I'/ujg = 0.1 and Us = c/10. The spatial distance x from the shock at x = is measured in CR Larmor radii. The vertical axis is scaled 
such that /jj" = 1 at the shock. 



ferent as seen in Figure 2(c) and in figure 3. Diffusion theory 
can still be applied to CR at a perpendicular shock by de- 
constructing the CR trajectory into difTusive motion of the 
gyrocentre and gyration about the gyrocentre. The gyratory 
part of the trajectory acts as a spatial convolution on the 
motion of the gyrocentre. This spreads the CR trajectory 
over a distance of one Larmor radius and smooths the CR 
density at the shock. A discontinuous jump in CR density 
is not possible at the shock on a scale less than a Larmor 
radius. Correspondingly, figures 2(c) and 3 exhibit a CR 
density gradient both downstream and upstream of a per- 
pendicular shock. Consequently, the CR density precisely at 
the shock is smaller than the CR density far downstream. 
In this case, the rate at which CR cross the shock is reduced 
and the spectrum at quasi-perpendicular shocks is steeper 
than p"* in line with equations (7) and (8). According to 
equations (9) and (10) of Jokipii (1987) the diffusive approx- 
imation is only valid at a quasi-perpendicular shock when 
i^/iOg > Us/c. This is consistent with the results shown in our 
figure 1 where we see that deviations from a spectrum 
occur where this condition is not met. 

The profile for 6 — 72° is intermediate between the 
quasi-perpendicular and the quasi-parallel cases. The steep- 



ening and flattening effects are both present but cancel out 
to leave a CR density approximately equal to that far down- 
stream and a CR spectrum proportional to p~*. 



5 OBSERVATIONS 

(i) Galactic cosmic rays 

CR arriving at the earth with energies up to a few PeV 
conform to a remarkably straight power law spectrum p ' 
with 7 equal to 4.6-4.7. The more rapid escape from the 
galaxy of high energy CR steepens the spectrum, but this 
cannot be solely responsible for steepening the spectrum 
from an index at source of 7 = 4.0. The underlying CR 
source spectrum probably has an index between about 4.2 
and 4.3 (Gaisser, Protheroe & Stanev 1998, Hillas 2005). 
Figure 1 indicates that the oblique-shock effects discussed 
here can produce a steepened spectrum (7 > 4) at high 
shock velocities, but only if the shock has a tendency towards 
a quasi-perpendicular rather than quasi-parallel configura- 
tion. As discussed in section 4, the spectral index averaged 
over random magnetic field orientations is approximately 
the same as for the low velocity limit: 7 = 4. 

The shock might be preferentially quasi-perpendicular 
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Figure 3. Spatial profiles of the CR density /g within one CR Larmor radius either side of the shock, u/uig = 0.1 and Us = c/10. 



in the event of strong magnetic field amplification in the 
precursor due to CR streaming (Bell 2004, 2005). Field am- 
plification is driven by jcfl x B where ]cr is the electric 
current carried by streaming CR and B is the local mag- 
netic field. The plasma is set in motion perpendicular to 
jcR and therefore predominantly perpendicular to the shock 
normal since jcR is directed away from the shock. These up- 
stream fluid motions stretch the magnetic field perpendic- 
ular to the shock normal producing a quasi-perpendicular 
shock orientation as evident in the sample field lines shown 
in figure 4 of Bell (2005). Hence strong field amplification can 
be expected to lead preferentially to a quasi-perpendicular 
shock geometry and the steepening of the spectral index. 
This is not contradicted by polarimetric observations which 
suggest a predominantly radial magnetic field in SNR. Zi- 
rakashvili & Ptuskin (2008) point out that this field orienta- 
tion probably develops downstream of the shock due to dif- 
ferential fluid motions after the plasma has passed through 
the shock. A similar effect is seen in simulations by Giacalone 
& Jokipii (2007). Alternatively, the shock might be prefer- 
entially quasi-perpendicular if the supernova expands into a 
pre-existing circumstellar wind supporting a magnetic field 
in the form of a Parker spiral. This would be unsurprising 
since SN are often preceded by strong mass loss. 

Another possible source of spectral steepening is non- 
linear feedback onto the shock hydrodynamics (Drury & 
Volk i98i, Drury 1983). Non-linear effects produce a con- 
cave spectrum which is steepened at low energy but fiatter 
than p"* at high energy (Bell i987, Falle & Giddings i987). 
Allen et al (2008) find evidence for concavity in the syn- 
chrotron spectrum of SNi006. The galactic CR spectrum is 
not noticeably concave, but a compensating convexity due 
to momentum-dependent escape from the galaxy could pro- 
duce the measured spectrum. The CR spectrum could be the 
result of a complicated mixture of non-linear, oblique-shock 
and escape effects. 



(ii) Young supernova remnants 

Further evidence for deviation from a CR source 
spectrum can be found in synchrotron spectra. SNR ra- 
dio spectral indices can deviate considerably from a = 0.5 
(7 — 4). Since synchrotron loss times at radio wavelengths 
are longer than the dynamical time, the measured spectral 
index gives a proper indication of the cosmic ray source spec- 
trum. Observations show that young pre-Sedov SNR con- 
sistently exhibit radio spectral indices significantly steeper 
than a = 0.5. Either magnetic field amplification or expan- 
sion into a Parker spiral may produce a geometry favouring 
quasi-perpendicular shocks and spectral steepening. But in 
either case, a high shock velocity is needed. If this is the 
explanation of the steeper than expected galactic CR spec- 
trum, then it implies that galactic CR are accelerated by 
young SNR rather than by SNR in the Sedov phase. 

In Figure 4 we plot a sample of SNR (see Table [T]) and 
their spectral indices as a function of shock velocity. It is re- 
lated to the analysis of Glushak (i985) who plotted the spec- 
tral indices of young shell SNR against SNR age to show that 
the spectra fiattened with time. Here we plot the spectral 
index against the mean expansion velocity Vsh calculated as 
the mean radius divided by SNR age, and extend and update 
the sample of SNR. We take the mean rather than the mo- 
mentary expansion velocity since the overall spectrum is an 
addition of cosmic rays that were accelerated throughout the 
entire age of the remnant. Depending on the number of cos- 
mic rays that enroll in the acceleration process early versus 
late, additional differences in the spectral index can result. 
The average time since acceleration would be different for 
core-collapse SNR and Type la SNe. In Type la SNe most 
CR are accelerated at late times as the shock expands into 
a uniform medium, whereas CR acceleration is more evenly 
spread in the case of core-collapse SNR which expand into 
a r~^ density profile of a circumstellar wind (Schure et al 
2010). 
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The Galactic remnants in figure 4 were taken from tfic 
updated online catalogue by Green (2009). They were se- 
lected if they were classified as having a clear shell morphol- 
ogy (i.e. no contamination in radio from a central source), 
had a measured spectral index and a readily available dis- 
tance estimate. Only sources were selected that are esti- 
mated to be younger than 10,000 yr, and not clearly in- 
teracting or overlapping with other radio sources. 

Because of a lack of data for galactic SNR with ages 
less than 200 years, we turn to the extragalactic remnants 
of recent SN for data on very young SNR. We added data 
for SN1987A (Manchester et al 2002, Manchester et al 2005) 
and SN1993,J (Wciler et al 2007), and a sample of Type 
Ib/c supcrnovac taken from Chevalier & Fransson (2006). 
Two data points are plotted for SN1987A corresponding to 
observations in 2002 and 2005 since its expansion velocity 
slowed significantly between these years. 

For our purpose, very young SNR from Type Ib/c SNe 
are especially interesting because of their high expansion ve- 
locities. They do not necessarily reflect a statistical sample. 
Their expansion velocities are based on models of the radius 
of the remnant at the time of the peak in synchrotron self 
absorption, as described in (Chevalier 1998). Each of the 
Type Ib/c's has a steep spectral index except for SN1998bw 
which was the only one to be associated with a GRB. Cheva- 
lier & Fransson (2006) argue that the fiattening with time 
may be a result of the nonlinear feedback of the cosmic rays 
in the plasma, which weakens with time. Nonlinear mod- 
els by Ellison & Reynolds (1991) can partly account for the 
steepening, up to a w 0.8, if the cosmic ray pressure con- 
stitutes a dominant fraction. The total steepening may be 
a cumulative effect of non-linear and oblique-shock steepen- 
ing. 

Except in the cases of the historical supernova rem- 
nants, the age of the galactic remnants is a highly uncertain 
parameter. Often it is based on the distance measurement, 
giving the size of the remnant, and from the temperature 
in X-rays a shock velocity is derived. Subsequently a Sedov 
model is used to arrive at the age estimate. However, in ad- 
dition to uncertainties in the size, if any energy is lost in 
e.g. cosmic rays the shock velocity may be severely underes- 
timated, as was shown might be the case in RCW 86 (Helder 
et al 2009). The data points in Figure 4 should therefore 
be considered with large error bars in Vsh and are likely to 
underestimate the shock velocity on average. The ones for 
which this is known arc indicated with an arrow, however 
this does not mean that other data points are necessarily 
more certain. Occasionally the shock velocity can possibly 
be overestimated when the age measurement is based on the 
age of the associated pulsar (i.e. G292.2-0.5). The age of the 
historical remnants is obviously more certain. However, their 
distance estimate may still introduce a large uncertainty in 
the size and shock velocity. Nevertheless, the data show a 
clear trend. As a guide, we plot on figure 4 a trend-line: 

a = 0.7-h0.31ogio(wsh/10'*km s"^) 

where Vsh is the radius divided by age for SNR. Although the 
data are uncertain at high velocity, the rise in the spectral 
index levels off at about a = 1. A levelling off also occurs 
at a ~ 0.5 for older SNR with lower expansion velocities. A 
better fit to the overall data may be an S-shaped curve with 
a stronger increase of spectral index at velocities around 



lO'^km s ^ . An appropriate trend- line for the historical SNR, 
for which the velocity estimates are more reliable, is 

a = 0.7 -I- 0.8 logio(wBh/10*km s"^) 

as plotted in figure 4. 

The two data points for SN1987A show that the varia- 
tion of spectral index with shock velocity is exhibited within 
the lifetime of a single SNR. The radio spectrum of Cas- 
siopeia A has similarly been observed to flatten with time as 
already argued in the seventies (Baars et al 1977). A flatten- 
ing of the Gas A radio spectrum from a = 0.79 to a = 0.77 
was found over a period of 15 years (1965-1980). Further 
flattening was conflrmed by O'Sullivan & Green (1990), al- 
though they found it proceeding at a slower rate. 

(in) Spectral index variations within SNR 

Further signatures of spectral index variation with 
shock obliquity may be visible in SNR expanding at high 
velocity into a uniform magnetic field. If magnetic field am- 
plification is weak, the field orientation at the shock is that 
of the pre-existing interstellar field. The electron spectrum 
should be flatter where the shock is quasi-parallel, increas- 
ing the number of CR electrons accelerated to high energy, 
and producing stronger X-ray synchrotron emission at the 
poles. 

Detailed, resolved, spectra arc available for the remnant 
of SN1006. This image shows a distinct bipolar symmetry, 
which is thought to reflect the directionality of the ambi- 
ent magnetic field. Most likely the bright caps are the parts 
where the shock velocity is parallelly aligned with the back- 
ground magnetic field (Rothcnflug et al 2004, Bocchino et 
al 2011). In an analysis of the combined radio- with the 
X-ray spectrum Allen et al (2008) found evidence for hard- 
ening of the spectrum in this area, both for models which 
included curvature in the spectrum and those without. The 
fractional hardening is modest but significant. The magnetic 
field in SN 1006 is generally believed to be not much am- 
plified and therefore plausibly reflects the original field ge- 
ometry. The spectral steepening towards the fainter parts, 
i.e. higher obliquity angle between the shock normal and 
magnetic field lines, is consistent with our findings. 

Detailed data are also available for Tycho's SNR, as was 
presented in Katz-Stone (2000). This remnant shows a much 
more patchy profile in radio. Spectra were extracted from 
various parts over the remnant. No clear trend with position 
was found. However, it was argued that the brightest parts 
have a flatter spectral index than the steeper parts, which 
is a similar trend to that in SN1006. 

A spatially resolved study of the radio spectrum of Cas- 
siopeia A (Wright et al 1999) shows variation of the spec- 
tral index with position in the remnant. Generally, a steeper 
power law is found in the outermost regions, whereas in the 
inner bright ring the spectral index is flatter (although still 
a > 0.75). The spectra they found do not show signs of spec- 
tral curvature and seem to be coherent over larger scales 
than the diffusion scale. This indicates a common mecha- 
nism, such as e.g. a preferentially perpendicular field at the 
forward shock might be able to account for. The magnetic 
field in Cas A is strongly amplified and thus a preferentially 
perpendicular field may be expected. 

(iv) Older supernova remnants 
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The radio spectral indices of old SNR in the self-similar 
Sedov phase are usually around 0.5, but many have flatter 
spectra which may be due to additional second order Fermi 
acceleration in MHD turbulence (Ostrowski 1999), or possi- 
bly due to compression ratios greater than four at radiative 
shocks. Because the shock velocities are much lower than in 
young SNR, the oblique-shock effects considered in this pa- 
per are unlikely to contribute to the flattening of the spectral 
index. 



Community's Seventh Framework Programme (FP7/2007- 
2013) / ERC grant agreement no. 247039 and from grant 
number ST/H001948/1 made by the UK Science Technology 
and Facilities Council. We especially thank John Kirk of 
MPIK Heidelberg for very helpful discussions. 



6 CONCLUSIONS 

We have shown that a wide range of spectral indices are pos- 
sible for shock velocities of around c/30 or higher. The CR 
spectrum is steepened for quasi-perpendicular shocks and 
flattened for shocks closer to parallel geometry. This may 
present a natural explanation of why few observed energetic 
particle spectra conform to the non-relativistic test-particle 

7 = 4 spectrum, although non-linear effects and second or- 
der Fermi acceleration probably also play a role. For a shock 
propagating into a randomly orientated upstream magnetic 
field the steepening and flattening of the spectrum averages 
out to produce 7 = 4. Observed radio spectra show that the 
spectrum is steeper at high shock velocity which suggests 
that high velocity shocks tend to be quasi-perpendicular 
rather than parallel. A predominantly quasi-perpendicular 
configuration may arise if magnetic field amplification is 
strong as expected at high shock velocity. Alternatively, 
a quasi-perpendicular geometry may result from expansion 
into a circumstellar wind supporting a Parker spiral. 

The galactic CR energy spectrum and the radio spec- 
tra of young SNR require a CR spectrum at source which is 
steeper than p"**. The galactic CR spectral index may be ex- 
plained if CR are accelerated by shocks such as those found 
in the historical SNR. These SNR have a relatively large 
shock area, still expand at a relatively high velocity into a 
possibly dense circumstellar medium, and produce CR with 
a suitable spectral index if the SNR shock is predominantly 
quasi-perpendicular. The spectral index is sensitive to the 
CR scattering rate i' which is presently not well determined 
by observation or theory. Observations (figure 4) suggest a 
fairly clear relation between SNR shock velocity and spectral 
index. 

Our model depends on the presumed simplification that 
the magnetic field structure can be separated into a small 
scale disordered component that scatters CR through small 
angles, and a large scale component which is uniform over a 
CR Larmor radius. However, the essential reason for spectral 
steepening by quasi-perpendicular shocks is that the mag- 
netic field sweeps CR through the shock making it difficult 
for CR to return to the shock for further acceleration. This 
effect can be expected to be independent of the detailed field 
structure. Another realisation of a related effect would be 
spectral steepening due to field line wandering (Kirk, Duffy 
& Gallant 1996). 
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Figure 4. Synchrotron spectral index versus stiock velocity for galactic (black diamonds), historic (blue filled circles), and extragalactic 
SNR and SNe (red triangles), as given in Table 1. The arrows indicate known issues with distance and / or age estimates and the likely 
limit derived from that in shock velocity. Note that there remain large errors in Dgji also for the other data points. The dashed line shows 
the trend in spectral index that is biased towards the historical galactic remnants a = 0.7 + 0.81ogiQ(tisij/10''km s~^). The dotted line 
shows a more average trend: a = 0.7 + 0.3 logiQ(Dsh/10'*km s^^). 
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Table 1. Overview of the supernova remnant models. Data are 
taken from the updated Galactic SNR catalogue by Green (2009) 
unless other source specified. Other sources: ''Gomez & Rodriguez 
(2009), ^Sankrit et al (2005), Chiotellis et al (2011), '^Zhou et 
al (2009), =Tian & Leahy (2008), ■''Chen et al (2004), SHarrus 
& Slane (1999), ''Hwang et al (2000), 'Haslam et al (1980), 
JGreiner et al (1994), ^Crawford et al (2001), 'Gaensler et al 
(1998), '"Gotthelf et al (1997), "Rakowski et al (2001), ^Slane et 
al, "JKinugasa et al (1998), ''Manchester et al (2002), ^Manchester 
et al (2005), *Weiler et al (2007), "Chevalier & Fransson (2006) 



